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Abstract  -  The  simulation  of  the  propagation  of 
electrical  activity  in  a  realistic-geometry  computer 
model  of  the  ventricles  of  the  human  heart  using  the 
governing  reaction-diffusion  equation  is  described. 
Each  model  point  is  represented  by  the  phase  1  Luo- 
Rudy  membrane  model,  appropriately  modified  to 
represent  human  action  potentials.  A  separate 
longer-duration  action  potential  waveform  was  used 
for  the  M  cells  found  in  the  ventricular  mid-wall. 
Cardiac  fiber  rotation  across  the  ventricular  wall 
was  implemented  via  an  analytic  equation,  resulting 
in  a  spatially-varying  anisotropic  conductivity  tensor 
and  consequently  anisotropic  propagation.  Since  the 
model  comprises  approximately  12  million  points, 
parallel  processing  was  used  to  cut  down  on 
simulation  time.  The  model  generated  acceptably- 
normai  electrocardiograms,  vectorcardiograms  and 
body  surface  potential  maps  on  the  surface  of  a 
numerical  human  torso  model.  Interestingly,  it  was 
found  that  the  intrinsic  difference  in  action  potential 
duration  between  M  cells  and  other  myocardial  cells 
was  greatly  diminished  due  to  electrotonic  coupling. 
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i.  Introduction 

Computer  models  have  long  been  used  for  the 
simulation  of  electrical  activity  in  the  heart.  Many  of 
these  propagation  models  were  of  the  cellular 
automaton  type,  simulating  the  electrical  activity  of  the 
model  cells  with  a  ruled-based  algorithm,  e.g.  Lorange 
and  Gulrajani  [1],  An  exception  was  the  model  of  Leon 
and  Horacek  [2]  that  used  subthreshold  electrotonic 
conduction  to  bring  the  cells  to  threshold,  with 
subsequent  action  potential  generation  and  refractory 
period  determination  being  ruled-based.  Only 
Huiskamp  [3]  recently  described  a  ventricular  heart 
model  in  which  both  subthreshold  and  suprathreshold 
activity  were  determined  by  modified  Beuler-Reuter 
membrane  equations  [4].  We  present  here  a  similar 
heart  model  in  which  each  model  cell  is  represented  by 
a  modified  version  of  the  more  accurate  phase  1  Luo- 
Rudy  [5]  membrane  equation.  Since  the  spatial 
resolution  is  greater  than  that  of  Huiskamp ’s  model 
(0.25  mm  as  opposed  to  0.6  mm),  simulations  with  the 
much  larger  number  of  points  (approximately  12  million 
as  opposed  to  800,000  in  Huiskamp’s  model) 


necessitated  parallel  processing  on  a  multiprocessor 
computer. 


n.  Methodology 


The  resolution  of  the  problem  was  done  in  two 
steps:  the  propagation  of  the  electrical  activity  in  the 
heart  was  simulated  first,  then  the  determined 
transmembrane  potentials  were  used  to  calculate  the 
torso  surface  potentials. 

A.  Simulation  of  Propagation 

The  propagation  of  electrical  activity  was  simulated 
in  the  ventricles  of  the  human  heart  model  developed  by 
Lorange  and  Gulrajani  [1],  The  ventricular  myocardium 
can  be  represented  by  two  continuous  domains, 
intracellular  and  interstitial,  characterized  by  the 
equations, 
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respectively.  Here  4>i  and  <()c  are  the  intracellular  and 
interstitial  potentials,  respectively,  Gj  and  Ge  the 

intracellular  and  interstitial  conductivity  tensors,  (3  the 
surface  to  volume  ratio  of  the  cardiac  cells,  Istim  the 
intracellular  stimulation  current,  and  Im  the 
transmembrane  current  coupling  the  two  domains.  Im 
passes  from  the  intracellular  to  the  interstitial  space  and 
is  the  sum  of  the  ionic  and  capacitive  currents: 
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where  Cm  is  the  specific  membrane  capacitance  (Cm  =  1 
pF/cm2)  and  Vm  is  the  transmembrane  potential  (V,„=  4>; 
-  cj)c).  The  combination  of  equations  (1),  (2)  and  (3)  and 
the  approximation  of  equal  anisotropy  in  the 
intracellular  and  interstitial  spaces  (G;  =  cGe )  results  in 
the  reaction-diffusion  equation: 
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As  mentioned,  the  dynamics  of  the  ionic  currents 
were  based  on  the  Luo-Rudy  [5]  membrane  model.  The 
action  potential  duration  was,  however,  modified  to 
match  that  of  human  ventricular  cells.  Two  different 
action  potential  durations  were  used,  a  longer  one  for 
the  so-called  M  cells  [6]  situated  in  the  mid-ventricular 
walls  and  a  shorter  one  for  endocardial  and  epicardial 
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cells.  These  changes  were  realized  by  decreasing  the 
time  constants  id  and  if,  responsible  for  the  slow  inward 
current  Is;  by  a  factor  of  1.6  for  the  M  cells,  and  by  a 
factor  of  2.0  for  the  endocardial  and  epicardial  cells. 
Furthermore,  the  potassium  conductance  GK1, 
responsible  for  the  time-independent  potassium  current 
IKi,  was  diminished  for  the  M  cells.  These  changes 
produce  a  transmural  gradient  of  transmembrane 
potential  during  repolarization  that  contributes  to  the 
upright  T  waves  noted  in  the  normal  electrocardiogram 
(ECG). 

The  propagation  of  electrical  activity  starts  at  the 
transition  points  between  the  Purkinje  fibers  and  the 
myocardium;  these  points  are  stimulated  at 
predetermined  times.  The  propagation  is  more  rapid 
along  the  cardiac  fibers  since  myocardial  anisotropy  is 
included  in  the  model.  The  fiber  orientations  are  defined 
analytically  with  a  modification  of  the  Beyar  and 
Sideman  [7]  equation.  The  longitudinal  and  transverse 
interstitial  conductivities  along  and  across  fiber 
directions  are  2.22  mS/cm  and  1.33  mS/cm  [8]. 
Intracellular  conductivities  are  half  these  values  as  2, 
was  taken  as  0.5. 

Finite  differences  were  used  to  solve  equation  (4). 
The  temporal  component  was  solved  with  the  forward 
Euler  method.  The  method  of  Victorri  et  al.  [9]  and  the 
use  of  table  look  ups  permitted  rapid  determination  of 
the  Luo-Rudy  ionic  currents.  The  formulation  of 
Salaheen  and  Ng  [10]  was  used  for  the  calculation  of  the 
anisotropic  diffusion  current  given  by  the  first  term  on 
the  right  in  equation  (4).  A  time  step  of  25  ps  was  used 
during  the  action  potential  upstroke;  this  was  increased 
to  50  ps  during  the  plateau.  The  crossover  from  one 
time  step  to  the  other  occurred  when  all  model  points 
terminated  their  upstrokes.  The  value  of  P  was  adjusted 
to  control  the  total  excitation  time.  Table  I  presents  the 
parameters  used  to  solve  equation  (4). 


TABLE  I 

Simulation  Parameters 


Parameter 

Symbol 

Value 

Stimulus  current 

Istim 

200  |aA  cm'2 

Spatial  step 

Ax 

250  nm 

Membrane  capacitance 

Cm 

1  |rF  cm'2 

Surface  to  volume  ratio 

P 

333  cm'1 

Equal-anisotropy  factor 

0.5 

Longitudinal  interstitial  conductivity 

gel 

2.22  mS  cm'1 

Transverse  interstitial  conductivity 

get 

1.33  mS  cm'1 

Small  time  step 

At, 

25  |us 

Large  time  step 

At2 

50  |us 

The  propagation  of  electrical  activity  was  simulated 
on  a  Silicon  Graphics  Origin  2000  parallel  computer 
consisting  of  64  R12000  processors  operating  at  400 
MHz.  Running  the  simulation  on  the  8  processors 
accessible  to  us  took  20  hours.  The  use  of  finite 
differences  allowed  an  easy  separation  of  the  problem 
into  portions  that  could  be  handled  simultaneously  on 
different  processors.  Total  simulation  time  depends 


mainly  on  the  time  spent  to  calculate  the  nodal  ionic 
currents,  and  using  the  full  bank  of  64  processors  would 
likely  cut  simulation  time  by  very  nearly  a  factor  of  8. 

B.  Torso  Surface  Potentials 

The  spatial  gradient  of  the  transmembrane  potential 
distribution  determined  from  equation  (4)  was  used  to 
calculate  elemental  current  dipoles  at  each  model  point. 
This  assumes  the  approximation  of  an  isotropic 
myocardium,  and  results  in  dipoles  that  are  everywhere 
normal  to  the  excitation  wavefront.  The  elemental 
dipoles  were  combined  into  58  regional  dipoles  that,  in 
conjunction  with  a  human  torso  model,  were  then  used 
to  calculate  the  torso  surface  potentials  using  standard 
integral  equation  approaches  [11]. 

hi.  Results 


Following  initial  adjustment,  the  heart  model 
generated  an  acceptably-normal  ECG  (Fig.  1).  Also 
verified  to  be  normal  were  the  vectorcardiogram  and  the 
body  surface  potential  map.  Note  the  largely  upright  T 
waves  in  the  ECG  of  Fig.  1,  and  as  mentioned  this  was 
realized  largely  by  the  presence  of  the  M  cells  in  the 
ventricular  walls,  with  their  longer  action  potential 
durations. 
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Fig.  1:  The  simulated  normal  ECG.  Depolarization 
(QRS)  and  repolarization  (T  wave)  are  both  shown. 


Fig.  2:  Action  potential  waveforms  for  A)  endocardial 
or  epicardial  cells  and  B)  M  cells. 


TABLE  II 

Action  Potential  Duration  Across  the  Left 
Ventricular  Wall 


Endocardial  cell 

M  cell 

Epicardial  cell 

274  ms 

310  ms 

255  ms 

Fig.  2  shows  the  intrinsic  action  potential  durations 
of  the  M  cells  as  well  as  those  of  the  endocardial  or 
epicardial  cells.  There  is  an  intrinsic  difference  in 
action  potential  duration  of  approximately  80  ms. 
During  simulation,  however,  this  intrinsic  difference 
was  greatly  diminished  by  the  electrotonic  coupling 
introduced  between  adjacent  cells  by  equation  (4). 
Table  II  gives  the  action  potential  durations  measured  at 
one  site  across  the  left  ventricular  wall  during  normal 
excitation.  It  is  seen  that  the  intrinsic  difference  in 
action  potential  durations  drops  from  80  ms  to 
approximately  35-55  ms. 

iv.  Discussion 

We  have  demonstrated  the  feasibility  of  simulating 
normal  activation  of  the  ventricles  using  the  reaction- 
diffusion  equation  (4),  together  with  an  accurate 
membrane  model  for  the  ventricular  cells.  The 
simulated  ECG  of  Fig.  1  is  close  to  a  normal  ECG. 
Even  more  interesting  is  the  demonstration  of  the 
reduction  in  intrinsic  action  potential  durations  across 
the  ventricular  wall  due  to  electrotonic  coupling. 
Anyukhovsky  et  al.  [12]  noted  that  the  difference  in 
action  potential  durations  between  M  cells  and 
endocardial/epicardial  cells  was  greater  in  vitro  than  in 
situ,  and  suggested  that  this  was  due  to  electrotonic 
coupling.  Verification  of  this  electrotonic  coupling 
hypothesis  is  difficult  experimentally,  but  as  evidenced 
by  the  above  results,  is  easy  to  accomplish  via  computer 
simulation. 


In  order  to  cut  down  on  computation  time,  we  used 
the  simple  criterion  of  switching  from  the  25  ps  to  the 
50  ps  time  step  once  the  upstroke  terminated  at  all 
model  points.  This  criterion  is  easily  implemented  for 
the  simulation  of  a  single  normal  beat.  When  simulating 
arrhythmias,  however,  multiple  disparate  wavefronts 
from  consecutive  beats  may  coexist  at  any  given  time 
instant,  and  a  more  sophisticated  space -time  wavefront 
tracking  algorithm  is  needed  to  determine  where  and 
when  the  time  step  needs  to  be  switched.  One  such 
adaptive  tracking  algorithm  was  proposed  recently  by 
Cherry  et  al.  [13],  and  the  incorporation  of  such  an 
algorithm  into  the  model  is  essential  in  order  to  render  it 
truly  versatile. 

v.  Conclusions 

The  use  of  the  more  physiologically-correct 
reaction-diffusion  equation  (4)  and  the  Luo-Rudy 
membrane  model  [5],  together  with  two  types  of  action 
potential  waveforms,  allows  for  more  realistic  ECG 
simulations.  This  kind  of  large-scale  simulation  is  only 
feasible  with  parallel  processing  rather  than  with 
conventional  serial  computing.  We  have  demonstrated 
the  computation  of  the  normal  ECG.  The  next  step  is  to 
implement  a  more  sophisticated  wavefront  tracking 
algorithm  into  the  model  in  order  to  be  able  to  simulate 
cardiac  arrhythmias. 
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